Dynamics of Freely Cooling Granular Gases 
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We study dynamics of freely cooling granular gases in two-dimensions using large-scale molecular 
dynamics simulations. We find that for dilute systems the typical kinetic energy decays algebraically 
with time, E(t) ~ t _1 , in the long time limit. Asymptotically, velocity statistics are characterized 
by a universal Gaussian distribution, in contrast with the exponential high-energy tails character- 
izing the early homogeneous regime. We show that in the late clustering regime particles move 
coherently as typical local velocity fluctuations, Av, are small compared with the typical velocity, 
Av/v ~ t -1//4 . Furthermore, locally averaged shear modes dominate over acoustic modes. The small 
thermal velocity fluctuations suggest that the system can be heuristically described by Burgers-like 
equations. 
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Freely evolving granular media, i.e., ensembles of hard 
sphere particles undergoing dissipative inelastic colli- 
sions, exhibit interesting collective phenomena includ- 
ing clustering, vortices, shocks, and multiple dynamical 
regimes. This dissipative nonequilibrium gas system can 
be described by kinetic theory in the early homogeneous 
phase where density fluctuations and velocity correla- 
tions are relatively small. However, quantitative charac- 
teristics of the late clustering regime such as the typical 
length and velocity scales, particle velocity distributions, 
as well the corresponding continuum theory remain open 
questions §||||gflg|j^Q||gg|||0|||| . 

Recent molecular dynamics simulations in one- 
dimension have shown that the asymptotic energy decay 
is universal and that clustering corresponds to the for- 
mation of shocks in the inviscid Burgers equation [fTS. 
In this paper we study long time asymptotic dynamics 
of freely evolving granular gases in two dimensions us- 
ing Molecular Dynamics (MD) simulations. Treating the 
particles as identical, totally undeformable hard disks, we 
developed an efficient event-driven algorithm pp| , pl| , p2| 
that allows us to probe the system well into the clus- 
tering regime. Our main result is that the asymptotic 
dynamics of dilute granular gases is universal as it is in- 
dependent of the degree of dissipation. 

In the simulations, N — 10 6 identical disks of ra- 
dius R — 0.15 and mass m = 1 are placed in a two- 
dimensional system of linear dimension L = 10 3 with 
periodic boundary conditions. The particle concentra- 
tion, c = 1, so the mass density (or volume fraction) is 
a = cttR 2 = 0.0707. Initially, particles are distributed 
randomly in space and their velocities are drawn from a 
Gaussian distribution with zero mean and unit variance. 
To ensure random initial conditions, the system is first 
evolved under elastic collisions only, with each particle 
undergoing 10 2 collisions on average. Then, time is reset 
to zero and each particle experiences up to an average of 
4 x 10 4 inelastic collisions. In such a collision the particle 
velocity, v, changes according to 
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(1 +r)(g • n)n, 



(1) 



with g = v — v' the relative velocity of the colliding par- 
ticles and n a unit vector connecting the centers of the 
two particles. In each collision, the normal component 
of the relative velocity is reduced by the restitution co- 
efficient, < r < 1, and the energy dissipation equals 
i(l — r 2 )(g • n) 2 . In the following, time is quoted in units 

°f ^° = 20cfl \Z m l which is proportional to the initial 
mean free time, where Eq is the average particle energy 
initially. 

Inelastic collapse, the formation of a cluster via an in- 
finite series of collisions occurring in a finite time, poses 
a difficulty for numerical simulations [Q . To properly re- 
solve such finite time singularities, we take collisions to 
be perfectly elastic (r — 1) when the relative velocity falls 
below a pre-specified threshold g • n < S. This scheme, 
which mimics actual granular particles where r —* 1 for 
sufficiently small relative velocities, can be applied as 
long as the cutoff velocity 6 is much smaller than the 
root mean square (rms) velocity, v lms = (w 2 ) 1 / 2 |p^| , p3| . 
In our simulations, i> rms > 10~ 3 , and 5 = 10 -5 . 

Given the initial conditions, the collision sequence is 
deterministic, and the system evolves freely (velocities 
are not controlled). In the absence of energy input, the 
system "cools" infinitely. The average particle energy, 
E(t) = \{v 2 ), as a function of time is shown in Fig. 1. 
First, we confirmed that for sufficiently small times, when 
the spatial distribution of particles remains roughly uni- 
form, the energy decay follows Haff 's cooling law 



E(t) =E (l 



(2) 

where g(r) is 



Here, (t /U)~ 1 = 3y/n{l - r 2 ]g(27?)/20, 
the pair correlation function JlQ JTTf] . 

However, this cooling law eventually breaks down due 
to the formation of dense clusters |,|,|jl^Jl|,|l|,|l|] . The 
deviation from this law occurs at a time scale that ulti- 
mately diverges when collisions become elastic. Interest- 
ingly, our simulations show that a universal decay law 
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E(t) ~ A(r)r 



(3) 



holds beyond this time scale. This implies that 
the typical velocity decays with time according to 
v ~ < -1 / 2 . The typical length scale explored by a par- 
ticle, C{t) ~ j Q dt'v{t') ~ t 1 ' 2 , remains small compared 
with the system size, C <C L, and thus, finite size effects 
were negligible throughout the simulations. The inset to 
Fig. f shows that the energy decay law is independent of 
the cutoff velocity as long as S < 1CP 4 (we also verified 
that numerical errors were irrelevant). 
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FIG. 1. The average particle energy E(t) as a function of 
time t for three different restitution coefficients, r — 0.1, 0.5, 
and 0.9. A line of slope —1 is shown for reference. The inset 
shows E(t) as a function of t for r = 0.1 and different cutoffs 
8 — 10 -4 , 10 -5 . Each line represents an average over four 
independent realizations. 

Overall, this behavior is consistent with the r- 
independent energy decay law, E ~ i~ 2 / 3 , found in one- 
dimension ||l2|| , although the exponent differs from the 
ID value of 2/3. This asymptotic law (||) has not been 
observed in previous numerical studies primarily due to 
insufficient temporal range p[j7|^|, pT|Jl3| , ^5[ . For exam- 
ple, power-law behavior was previously suggested by de- 
formable sphere molecular dynamics simulations in two- 
and three-dimensions jl3|. However, the value of the scal- 
ing exponent depended on parameters such as the system 
size and the restitution coefficient. 

The simulations show that as the volume fraction a 
was reduced, the dependence of the prefactor A(r) on the 
restitution coefficient r became weaker and weaker, sug- 
gesting completely universal behavior in the dilute limit, 
a — > 0. Further numerical studies with smaller volume 
fractions are needed to fully resolve this issue. If indeed 
the prefactor A is independent of the restitution coeffi- 
cient r, this implies a certain behavior of £ c , the time 
scale marking the transition from the homogeneous to 
the clustering regime. Matching the two asymptotic be- 
haviors E{t) ~ [(1 - r)t]- 2 for t < t c with E(t) ~ t' 1 



for t ^> t c shows that this crossover time scale diverges 
according to t c ~ (1 — r)~ 2 in the quasi-elastic limit, 
r -> 1. 

Next, we examined whether the entire velocity distri- 
bution, not merely the typical velocity scale is indepen- 
dent of the restitution coefficient. We studied P(v,t), 
the probability distribution function of the velocity mag- 
nitude v = |v| at time t. Figure 2 clearly shows that a 
universal scaling function underlies the velocity distribu- 
tion when t 3> t c (r) 



P(v,t) 
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This provides evidence that the asymptotic dynamics are 
universal. 

The simulations suggest that the corresponding scal- 
ing function is Gaussian, $(z) = 7r _1 exp(— z 2 ), (see 
Fig. 2). We therefore studied the kurtosis, particularly 
the quantity n — (w 4 )/(w 2 ) 2 — 2. Initially, k = 0, 
the expected value for a Gaussian distribution in two- 
dimensions. In the intermediate regime, k is distinctly 
nonzero and strongly depends on the value of the resti- 
tution coefficient. However, in the asymptotic regime 
(t t c ), k does not depend on restitution coefficient 
and is very close to zero < 0.1), a manifestation of 
the universal velocity statistics. Given this universality, 
we present data for a representative value of r = 0.1 in 
the rest of this paper. 
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FIG. 2. The scaling velocity distribution $(2:) versus the 
square of scaling variable z = v/v ms - The data corresponds 
to three different restitution coefficients r = 0.1, 0.5, and 0.9, 
at three different times, t = 5 x 10 2 (excluding r — 0.9), 
2 x 10 3 , 2 x 10 4 , all in the clustering regime. These eight 
distributions follow a universal scaling function, a Gaussian 
distribution. The inset shows the exponential high-energy tail 
of P(v) at t = 10 for r = 0.1, where the velocity was rescaled 
by the rms velocity. 

For completeness, we mention that in the intermediate 
homogeneous regime, the distribution has an exponential 
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high-energy tail, in agreement with kinetic theory studies 
P,p| Jl^jr7|fl8| , |l9| . The measured slope of 3.3 is consistent 
with Direct Simulation Monte Carlo (DSMC) @ and 
MD simulations (see the inset in Fig. 2). 
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FIG. 3. The adjusted kurtosis k as a function of time for 
three restitution coefficients r = 0.1, 0.5, and 0.9. 

Both the Gaussian and the exponential extremal ve- 
locity statistics are consistent with the following heuris- 
tic argument. The large- velocity tails are dominated 
by the fastest particles (v ~ 1) initially present in the 
system, still moving with their initial velocity at time 
t. These particles managed to avoid any collisions with 
other particles. The "survival" probability, S(v,t), for 
such particles decays exponentially with the volume they 
sweep till time t, S(v = oc exp(— const, x t). As- 
suming a stretched exponential large- velocity tail for the 
scaling function Q(z) ~ exp(— const, x |z| 7 ) as z — > oo 
and the typical velocity decay v ~ t _/3 as t — > oo, 
and then matching the dominant behavior of the fastest 
particles P(l,t) — P(l, 0)5(1, i) cx exp(— const, x t) 
with the expected scaling behavior P(l,t) oc $ (t^ oc 
exp (—const, x t' 37 ) yields the exponent relation /?7 = 1. 
In the homogeneous regime, (3 = 1 and thus 7=1, while 
in the clustering regime, (3=1/2 and thus, 7 = 2. 
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FIG. 4. The ratio of the thermal energy Eth to the kinetic 
energy Ek as a function of time for r = 0.1. A line of slope 
— 1/2 is also shown as a reference. 

To quantify collectives motions and the correspond- 
ing velocity fluctuations in the clustering regime, we also 



measured the local kinetic energy, Ek, and the local ther- 
mal energy, E t h- The local kinetic energy is defined 
by Ek = \{pv?) with the density p and the velocity 
u = (u x ,u y ) obtained by averaging the corresponding 
quantities over a small region of space. The local ther- 
mal energy Eth = h{( v x — u x ) 2 + (v y — u y ) 2 ) is obtained 
by subtracting the average local velocity u from the par- 
ticle velocity v . Space was divided into 128 x 128 small 
boxes. As shown in Fig. 4, initially, the thermal energy is 
large compared with the kinetic energy, indicating that 
particles essentially move independently and there are no 
collective motions. In contrast, in the asymptotic regime, 
t ^> t c (r = 0.1) ss 10 2 the thermal energy becomes much 
smaller than the kinetic energy, indicating that coherent 
motion of particles becomes dominant. The appearance 
of the collective motions was suggested by linear stabil- 
ity analysis of the hydrodynamic equatio ns an d ob- 
served in MD simulations |,|,|,@,^p,|lg|l|. 

Interestingly, we find that the ratio of thermal to ki- 
netic energy decays algebraically, E t h/Ek ~ t~ - 5 , in the 
long time limit. The thermal energy, E t h = ^((Af) 2 ) 
quantifies local velocity fluctuations, Aw, and using 
Eq. (1), 

E th ~(Avf~t-^. (5) 

In other words, local velocity fluctuations, Av ~ i~ 3 / 4 
are small compared with the typical velocity, v ~ t~ 1//2 , 
as Av/v ~ i -1 / 4 . Thus, at least two distinct velocity 
scales are needed to characterize particle velocities in the 
clustering regime. Moreover, the relatively small velocity 
fluctuations imply strong velocity correlations and a well- 
defined average local velocity. Thus, a hydrodynamic de- 
scription is plausible in the clustering regime. 
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FIG. 5. The ratio of the longitudinal energy Ei to the trans- 
verse energy E t as a function of time. 

To distinguish between acoustic and shear modes, we 
transformed the averaged velocity field into Fourier space 
and decomposed the velocity as u(k) = u;(k) + u t (k). 
Here u t (k) is the transverse velocity satisfying k • u t = 0, 
and u;(k) is the longitudinal velocity with k x u; = 0. 
The corresponding representation in physical space is 
u(x,y) = ui(x,y) + u t (x,y), with V • u t = 0, and 
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V x u; = 0. The averaged longitudinal energy, i.e., the 
energy of acoustic modes, is defined as the spatial av- 
erage Ei = -k(puf), and the transverse energy, i.e., the 
energy of shear modes, is defined as the space average of 
Et = h(pv,f). Initially, the system is isotropic and conse- 
quently, the longitudinal and the transverse energies are 
equal (see Figure 5). Since collective motions are negligi- 
ble initially, statistical fluctuations are isotropic as well. 
The ratio of longitudinal to transverse energies decreases 
steadily and eventually, it saturates at a value of roughly 
Ei/Et — > 0.4 for t > t c . Therefore, shear modes domi- 
nate over acoustic modes in the clustering regime. This 
behavior is consistent with the formation of dense, thin, 
elongated clusters, where particles move coherently, par- 
allel to the cluster orientation Interestingly, we did 
not observe obvious correlations between the cluster size 
and its velocity. 

Collective motion of granular media can be described 
by mass, momentum and energy balance equations 
P,[7| p^j25[ | . We have seen that the thermal energy, i.e., 
the temperature T, is negligible compared with the ki- 
netic energy, and therefore, we may expand the system 
around zero temperature. In particular, we can ignore 
the pressure term in the momentum equation due to the 
fact that P^T while keeping the viscosity term since 
v ~ T 1 / 2 M. The resulting governing equations are 



dtp + d a (pu a ) = 
dt(pu a ) + dp(pu a up) -- 

7Ta/3 = d a (pVUf}) + dp(pVU a ) 



0. 

= dp-Rap, 

- d 1 (pvu 1 )5 a p, 



(6) 
(7) 
(8) 



taken in the limit of vanishing viscosity, v — > 0. The 
above equations are similar to the two-dimensional Burg- 
ers equation supplemented by the continuity equation, 
used to model large-scale formation of matter in the 
universe ]26| , |27| ]. Here, the formation of shocks corre- 
sponds to dense, thin, string-like clusters where particles 
move parallel to the cluster orientation. At least qual- 
itatively, this is consistent with the above findings that 
shear modes dominate over acoustic modes. 

In conclusion, our simulations show that in the low 
density and large system limits, the kinetic energy of 
freely evolving granular media decays with time as i _1 
in the long time limit. We have also observed that the 
particle velocity distribution is Gaussian in this asymp- 
totic time regime. The above behavior is independent 
of the degree of dissipation. We have studied the origin 
of this universality and found that the asymptotic dy- 
namics are dominated by collective motions of particles 
or alternatively, shear modes. In the clustering regime, 
strong spatial velocity correlations develop, as local ve- 
locity fluctuations are much smaller compared with the 
typical velocity. We have also argued that the governing 
continuum equations are similar to the two-dimensional 
Burgers equation with infinitely small viscosity. 
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